Application of Control Charts in Pairs Trading: An Analysis of BIST30 Stock Indices

Project code, data results and data figures are below; discussion and conclusion are after these.

import pandas as pd
import plotly as plt
import numpy as np
import plotly.express as px
import os
import glob
from sklearn.linear_model import LinearRegression
import plotly.graph_objects as go

We get all data from excels and create a column for each stock.

path = os.getcwd()
csv_files = glob.glob(os.path.join(path, "*.csv"))

last_data = pd.DataFrame()

for f in csv_files:
    data = pd.read_csv(f)
    symbol = data.short_name.unique()
    new_data = pd.DataFrame()

    data['timestamp'] = pd.to_datetime(data['timestamp'])
    data['timestamp']  = data['timestamp'].dt.tz_convert('Europe/Istanbul')
    data.set_index('timestamp', inplace=True)
    data.index = data.index.tz_localize(None)
    #data.index = pd.to_datetime(data.index).strftime('%Y-%m-%d')

    for i in symbol:
        stock = data[(data['short_name'] == i)]
        stock_close = stock.groupby([stock.index]).agg('last')
        new_data[i]= stock_close['price']
    
    last_data = last_data.append(new_data)
last_data= last_data.dropna()

Calculation of correlation for first 3600 inputs, that is 1.5 year between 2018-January and 2019-July, and two of the most correlated stock pairs founded. These are GARAN-AKBNK and SAHOL-VKBNK.

correlations = last_data[:3600].corr()
correlations = correlations.replace(1, 0)
s = correlations.unstack()
so = s.sort_values(kind="quicksort")
so[-4:]

SAHOL  VAKBN    0.958566
VAKBN  SAHOL    0.958566
GARAN  AKBNK    0.964810
AKBNK  GARAN    0.964810
dtype: float64

fig = px.imshow(correlations, text_auto=True)
fig.update_layout(width=1000, height=1000)
fig.show()

This linreg function provide train data and call test_linear function to test linear regression.

def linreg(data, n1, n2, n3, symbol1, symbol2):
    df = pd.DataFrame()
    df[symbol1] = data[symbol1][n1:n3]
    df[symbol2] = data[symbol2][n1:n3]

    x_train = df[[symbol1]][:n2]
    y_train = df[symbol2][:n2]

    slr = LinearRegression()

    slr.fit(x_train, y_train)

    #prediction by train set
    y_expected_by_model= slr.predict(x_train)

    lindf = pd.DataFrame()
    lindf[symbol1+'_train'] = x_train[symbol1]
    lindf[symbol2+'_train'] = y_train
    lindf[symbol2+'_expected_by_model'] = y_expected_by_model
    lindf['Residuals'] = lindf[symbol2+'_train']-lindf[symbol2+'_expected_by_model']

    mean = lindf['Residuals'].mean()
    variance = lindf['Residuals'].var()
    UCL = mean + np.sqrt(variance)
    LCL = mean - np.sqrt(variance)

    test_data = test_linear(df, n1, n2, n3, slr, symbol1, symbol2)
    return lindf, UCL, LCL, mean, test_data

def test_linear(data_last, n1, n2, n3, linearregression, symbol1, symbol2):
    x_test = data_last[[symbol1]][n2:]
    y_test = data_last[symbol2][n2:]

    y_test_predict = linearregression.predict(x_test)

    test_df = pd.DataFrame()
    test_df[symbol1+'_test'] = x_test[symbol1]
    test_df[symbol2+'_test'] = y_test
    test_df[symbol2+'_predicted_by_model'] = y_test_predict
    test_df['Residuals'] = test_df[symbol2+'_test']-test_df[symbol2+'_predicted_by_model']

    test_df[symbol2] = data_last[symbol2][:]
    test_df[symbol1] = data_last[symbol1][:]

    return test_df

With symbols, we prepare data for buy and sell orders correctly by using two functions below. Our strategy is buying a stock if a residual is above UCL and selling if a residual is below LCL.

def order_and_profit(data, UCL, LCL, symbol1, symbol2):
    data[symbol2+'_buy'] = data[symbol2][data['Residuals']<=LCL]
    data[symbol1+'_buy'] = data[symbol1][data['Residuals']>=UCL]
    data = data.dropna(subset=[symbol2+'_buy', symbol1+'_buy'], thresh=1)
    
    data = data[(data[symbol2+'_buy'].notna() & data[symbol2+'_buy'].shift(1).isna()) | (data[symbol1+'_buy'].notna() & data[symbol1+'_buy'].shift(1).isna())]

    data[symbol1+'_sell'] = data[symbol1][data[symbol1+'_buy'].isna()]
    data[symbol2+'_sell'] = data[symbol2][data[symbol2+'_buy'].isna()]

    data_with_profit = profit_calculator(data, symbol1)
    data_with_profit = profit_calculator(data_with_profit, symbol2)

    data_with_profit = data_with_profit.drop(symbol1 + '_profit_sell', axis=1)
    data_with_profit = data_with_profit.drop(symbol2 + '_profit_sell', axis=1)

    return data
    

def profit_calculator(data, symbol):
    data[symbol + '_profit'] = (data[symbol + '_sell'] - data[symbol + '_buy'].shift(1))/data[symbol + '_buy'].shift(1)
    data[symbol + '_profit_sell'] = (data[symbol + '_sell'].shift(1) - data[symbol + '_buy'])/data[symbol + '_sell'].shift(1)

    data[symbol + '_profit'] = data[symbol + '_profit'].fillna(data[symbol + '_profit_sell'])
    data[symbol + '_profit'] += 1
    data.loc[data.index[0], symbol + '_profit'] = 100
    data[symbol + '_profit'] = data[symbol + '_profit'].cumprod()


    return data

We call functions for GARAN-AKBNK and test the trading strategy for the interval between 2019-July and 2020-July.

output = linreg(last_data, 0, 3600, 6000, 'GARAN', 'AKBNK')
linear_result = output[0]
LCL = output[2]
UCL = output[1]
mean = output[3]
test_df = output[4]

Linear Regression Model for GARAN-AKBNK with using first 3600 input as train data.

fig = go.Figure()

fig.add_trace(go.Scatter(x=linear_result['GARAN_train'], y=linear_result['AKBNK_train'],
                    mode='markers',
                    name='reality'))

fig.add_trace(go.Scatter(x=linear_result['GARAN_train'], y=linear_result['AKBNK_expected_by_model'],
                    mode='lines',
                    name='prediction'))
fig.show()

We use 1-sigma for limits and the graph is resulted like that with limits for train data. It is applied in linreg function.

fig = go.Figure()

fig.add_trace(go.Scatter(x=linear_result.index,
                        y=linear_result['Residuals'],
                        mode='markers',
                        name='reality'))

fig.add_hline(UCL)
fig.add_hline(LCL)
fig.add_hline(mean)

fig.show()

These are our orders and profits, GARAN creates 20% profit in 1 year and AKBKN creates a 1% loss in 1 year.

orders_and_profits = order_and_profit(test_df, UCL, LCL, 'GARAN', 'AKBNK')
orders_and_profits

	GARAN_test	AKBNK_test	AKBNK_predicted_by_model	Residuals	AKBNK	GARAN	AKBNK_buy	GARAN_buy	GARAN_sell	AKBNK_sell	GARAN_profit	GARAN_profit_sell	AKBNK_profit	AKBNK_profit_sell
timestamp
2019-07-01 10:00:00	8.7112	6.1421	6.434299	-0.292199	6.1421	8.7112	6.1421	NaN	8.7112	NaN	100.000000	NaN	100.000000	NaN
2019-09-18 17:00:00	8.4352	6.5104	6.240574	0.269826	6.5104	8.4352	NaN	8.4352	NaN	6.5104	103.168335	0.031683	105.996320	NaN
2019-12-13 12:00:00	9.8336	6.9216	7.222114	-0.300514	6.9216	9.8336	6.9216	NaN	9.8336	NaN	120.271735	NaN	99.301542	-0.06316

These graph is belong to test data, between 2019-July and 2020-July, and UCL and LCL are seen here.

fig = go.Figure()

fig.add_trace(go.Scatter(x=test_df.index, 
                        y=test_df['Residuals'],
                        mode='markers',
                        name='reality'))

fig.add_hline(UCL)
fig.add_hline(LCL)
fig.add_hline(mean)

fig.show()

Linear Regression Model for SAHOL-VAKBN with using first 3600 input as train data.

output1 = linreg(last_data, 0, 3600, 6000, 'SAHOL', 'VAKBN')
linear_result1 = output1[0]
LCL = output1[2]
UCL = output1[1]
mean = output1[3]
test_df1 = output1[4]

Linear Regression Model for SAHOL-VAKBN with using first 3600 input as train data.

fig = go.Figure()

fig.add_trace(go.Scatter(x=linear_result1['SAHOL_train'], y=linear_result1['VAKBN_train'],
                    mode='markers',
                    name='reality'))

fig.add_trace(go.Scatter(x=linear_result1['SAHOL_train'], y=linear_result1['VAKBN_expected_by_model'],
                    mode='lines',
                    name='prediction'))
fig.show()

We use 1-sigma for limits and the graph is resulted like that with limits for train data. It is applied in linreg function.

fig = go.Figure()

fig.add_trace(go.Scatter(x=linear_result1.index,
                        y=linear_result1['Residuals'],
                        mode='markers',
                        name='reality'))

fig.add_hline(UCL)
fig.add_hline(LCL)
fig.add_hline(mean)

fig.show()

These are our orders and profits, SAHOL creates 60% profit in 1 year and VAKBN creates a 10% loss in 1 year.

orders_and_profits = order_and_profit(test_df1, UCL, LCL, 'SAHOL', 'VAKBN')
orders_and_profits

	SAHOL_test	VAKBN_test	VAKBN_predicted_by_model	Residuals	VAKBN	SAHOL	VAKBN_buy	SAHOL_buy	SAHOL_sell	VAKBN_sell	SAHOL_profit	SAHOL_profit_sell	VAKBN_profit	VAKBN_profit_sell
timestamp
2019-06-28 16:00:00	7.2031	4.17	5.036464	-0.866464	4.17	7.2031	4.17	NaN	7.2031	NaN	100.000000	NaN	100.000000	NaN
2020-02-03 11:00:00	8.1708	6.59	6.221938	0.368062	6.59	8.1708	NaN	8.1708	NaN	6.59	86.565507	-0.134345	158.033573	NaN
2020-02-11 13:00:00	8.4316	6.18	6.541429	-0.361429	6.18	8.4316	6.18	NaN	8.4316	NaN	89.328551	NaN	167.865707	0.062215
2020-02-28 15:00:00	7.2199	5.43	5.057045	0.372955	5.43	7.2199	NaN	7.2199	NaN	5.43	102.165902	0.143709	147.493656	NaN
2020-03-05 10:00:00	8.2803	5.95	6.356080	-0.406080	5.95	8.2803	5.95	NaN	8.2803	NaN	117.171196	NaN	133.369033	-0.095764
2020-03-17 11:00:00	6.3111	4.30	3.943727	0.356273	4.30	6.3111	NaN	6.3111	NaN	4.30	145.036553	0.237817	96.384343	NaN
2020-04-29 15:00:00	7.1125	4.57	4.925476	-0.355476	4.57	7.1125	4.57	NaN	7.1125	NaN	163.453674	NaN	90.332303	-0.062791

These graph is belong to test data, between 2019-July and 2020-July, and UCL and LCL are seen here.

fig = go.Figure()

fig.add_trace(go.Scatter(x=test_df1.index, 
                        y=test_df1['Residuals'],
                        mode='markers',
                        name='reality'))

fig.add_hline(UCL)
fig.add_hline(LCL)
fig.add_hline(mean)

fig.show()

Moving Average Time Series Model for GARAN-AKBNK

def moving_average(data, symbol1, symbol2, w, n1, n3):
    df = pd.DataFrame()
    df[symbol1] = data[symbol1][n1:n3]
    df[symbol2] = data[symbol2][n1:n3]
    df[symbol1+'/'+symbol2] = df[symbol1]/df[symbol2]
    mean = df[symbol1+'/'+symbol2].mean()
    std = df[symbol1+'/'+symbol2].std()
    df['SMA30'] = df[symbol1+'/'+symbol2].rolling(w).mean()
    df['STD'] = df[symbol1+'/'+symbol2].rolling(w).std()
    df['UCL'] = df['SMA30']+df['STD']*3
    df['LCL'] = df['SMA30']-df['STD']*3
    df = df.dropna()

    return df

def order(data, symbol1, symbol2):
    data[symbol1+'_buy'] = data[symbol1][data[symbol1+'/'+symbol2]<=data['LCL']]
    data[symbol2+'_buy'] = data[symbol2][data[symbol1+'/'+symbol2]>=data['UCL']]

    data = data.dropna(subset=[symbol2+'_buy', symbol1+'_buy'], thresh=1)    
    data = data[(data[symbol2+'_buy'].notna() & data[symbol2+'_buy'].shift(1).isna()) | (data[symbol1+'_buy'].notna() & data[symbol1+'_buy'].shift(1).isna())]

    data[symbol1+'_sell'] = data[symbol1][data[symbol1+'_buy'].isna()]
    data[symbol2+'_sell'] = data[symbol2][data[symbol2+'_buy'].isna()]

    data_with_profit = profit_calculator(data, symbol1)
    data_with_profit = profit_calculator(data_with_profit, symbol2)

    data_with_profit = data_with_profit.drop(symbol1 + '_profit_sell', axis=1)
    data_with_profit = data_with_profit.drop(symbol2 + '_profit_sell', axis=1)

    return data_with_profit

GARAN/AKBNK parity is used for detection of mean and std for each 30 period. LCL is determined as mean minus 3std and UCL is determined as mean plus 3std. If data is below LCL, we buy GARAN and sell AKBNK; if data is above UCL, we buy AKBNK and sell GARAN.

GARAN creates 100% profit and AKBNK creates 40% loss in the same period as linear regression model.

output2 = moving_average(last_data, 'GARAN', 'AKBNK', 30, 3600, 6000)
profit_data = order(output2, 'GARAN', 'AKBNK')
profit_data

	GARAN	AKBNK	GARAN/AKBNK	SMA30	STD	UCL	LCL	GARAN_buy	AKBNK_buy	GARAN_sell	AKBNK_sell	GARAN_profit	AKBNK_profit
timestamp
2019-07-24 10:00:00	8.6578	6.4761	1.336885	1.324473	0.003669	1.335481	1.313465	NaN	6.4761	8.6578	NaN	100.000000	100.000000
2019-08-08 11:00:00	8.5420	6.4248	1.329536	1.369622	0.012827	1.408105	1.331140	8.5420	NaN	NaN	6.4248	101.337522	99.207857
2019-09-05 09:00:00	8.4975	6.3049	1.347761	1.338123	0.003198	1.347717	1.328528	NaN	6.3049	8.4975	NaN	100.809599	101.059280
2019-09-26 10:00:00	8.7112	6.7074	1.298745	1.317633	0.005781	1.334976	1.300289	8.7112	NaN	NaN	6.7074	98.274381	107.510827
2019-10-23 10:00:00	8.4885	6.0820	1.395676	1.345397	0.016745	1.395630	1.295163	NaN	6.0820	8.4885	NaN	95.762018	117.535169
2019-11-12 11:00:00	8.7557	6.3819	1.371958	1.384842	0.004157	1.397312	1.372372	8.7557	NaN	NaN	6.3819	92.747632	123.330762
2019-12-13 12:00:00	9.8336	6.9216	1.420712	1.400558	0.006013	1.418597	1.382520	NaN	6.9216	9.8336	NaN	104.165642	112.901013
2020-01-07 17:00:00	9.4417	6.6732	1.414868	1.428791	0.004276	1.441617	1.415964	9.4417	NaN	NaN	6.6732	108.316972	108.849260
2020-01-09 11:00:00	10.4660	7.1100	1.472011	1.433173	0.011882	1.468818	1.397528	NaN	7.1100	10.4660	NaN	120.067936	101.724439
2020-03-23 17:00:00	6.1637	4.6344	1.329989	1.396659	0.021521	1.461223	1.332094	6.1637	NaN	NaN	4.6344	169.424737	66.305449
2020-04-15 10:00:00	7.3930	4.9427	1.495741	1.471447	0.007674	1.494470	1.448425	NaN	4.9427	7.3930	NaN	203.215127	61.894528
2020-04-22 14:00:00	7.2772	4.8828	1.490374	1.506827	0.005263	1.522615	1.491039	7.2772	NaN	NaN	4.8828	206.398180	61.144436
2020-05-29 10:00:00	7.0456	4.8657	1.448014	1.424255	0.006855	1.444819	1.403690	NaN	4.8657	7.0456	NaN	199.829470	61.358569

fig = go.Figure()

fig.add_trace(go.Scatter(x=output2.index,
                        y=output2['GARAN/AKBNK'],
                        mode='markers',
                        name='reality'))

fig.add_trace(go.Scatter(x=output2.index,
                        y=output2['UCL'],
                        mode='lines',
                        name='UCL'))

fig.add_trace(go.Scatter(x=output2.index,
                        y=output2['LCL'],
                        mode='lines',
                        name='LCL'))

fig.show()

SAHOL creates 65% profit and AKBNK creates 64% loss in the same period as linear regression model.

output3 = moving_average(last_data, 'SAHOL', 'VAKBN', 30, 3600, 6000)
profit_data = order(output3, 'SAHOL', 'VAKBN')
profit_data

	SAHOL	VAKBN	SAHOL/VAKBN	SMA30	STD	UCL	LCL	SAHOL_buy	VAKBN_buy	SAHOL_sell	VAKBN_sell	SAHOL_profit	VAKBN_profit
timestamp
2019-07-04 11:00:00	7.5481	4.73	1.595793	1.654414	0.018318	1.709367	1.599461	7.5481	NaN	NaN	4.73	100.000000	100.000000
2019-10-01 10:00:00	8.0362	5.08	1.581929	1.531036	0.015088	1.576300	1.485772	NaN	5.08	8.0362	NaN	106.466528	92.600423
2019-11-13 16:00:00	7.2956	4.81	1.516757	1.547589	0.008232	1.572286	1.522892	7.2956	NaN	NaN	4.81	116.278269	87.678747
2020-02-04 16:00:00	8.5746	6.76	1.268432	1.243599	0.007891	1.267270	1.219927	NaN	6.76	8.5746	NaN	136.663145	52.133309
2020-05-08 11:00:00	7.0773	4.62	1.531883	1.563531	0.010302	1.594437	1.532626	7.0773	NaN	NaN	4.62	160.527317	35.629569
2020-05-14 11:00:00	7.2974	4.57	1.596805	1.565679	0.008573	1.591399	1.539959	NaN	4.57	7.2974	NaN	165.519625	36.015171

fig = go.Figure()

fig.add_trace(go.Scatter(x=output3.index, 
                        y=output3['SAHOL/VAKBN'],
                        mode='markers',
                        name='reality'))

fig.add_trace(go.Scatter(x=output3.index, 
                        y=output3['UCL'],
                        mode='lines',
                        name='UCL'))

fig.add_trace(go.Scatter(x=output3.index, 
                        y=output3['LCL'],
                        mode='lines',
                        name='LCL'))

fig.show()

DISCUSSION

In linear regression model, data between 2018 January and 2020 July is used. Data is divided as train data, from 2018 January to 2019 July, and test data, from 2019 July to 2020 July. Furthermore, the most correlated stocks are found for 2018 January/2019 July. Therefore, the model is not so successful in trading. In other words, it is assumed that we does not know anything about after 2019 July while creating trade simulation.

Simple trade simulation based on residuals between data predicted by linear regression model and real data. If residuals above UCL, we create a buy order for a stock and a sell order for another stock because they should be correlated for the model. We change orders residuals below LCL and try to catch cycles between stocks. However, the model does not catch that the correlation of stocks may decrease permanently; therefore, we may update model after every data comes true.

In moving average model, data between 2019 July and 2020 July is used because there is no need for train data. Also, stock1/stock2 parity is created. Mean and limits are updated for each data because time series should be adaptive for each new data. 30 is used for window and limits are created for 'mean+-std*3'. If the parity is above UCL or is below LCL, a buy order for a stock and a sell order for another stock are created until another out-of-limits situation.

CONCLUSION

Both methods do not create a great profit compared to normal cumulative return stocks for related dates. In stocks market, there are so external effects like political, sectoral and company related developments. Therefore, just looking correlation and relation between two data is wrong. However, our models can be better with different ordering buy-sell strategies and different control limit approaches like EWMA charts and so on.